package scu.maqiang.fes;

import static scu.maqiang.numeric.MVO.average;
import static scu.maqiang.numeric.MVO.matmul;
import static scu.maqiang.numeric.MVO.mul;
import static scu.maqiang.numeric.MVO.transpose;

import java.util.Arrays;

import scu.maqiang.mesh.*;
import scu.maqiang.mesh.Mesh2T6;
import scu.maqiang.numeric.CG;
import scu.maqiang.numeric.Constants;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.ParamCheck;
import scu.maqiang.numeric.SRMatrix;

public class FES2T62 extends FES {

	public FES2T62(Mesh2T6 mesh) {
		super(mesh, 2);
		eleStiff = FES2T62::elementElasticityStiff;
		eleStiff = FES2T62::elementElasticityMass;
		eleSource = FES2T62::elementVolumeForce;
		eleSource = FES2T62::elemenSurfaceForce;
	}

	
	/**
	 * ���㵯�����ⵥԪ�նȾ���
	 * 
	 * @param xy
	 *            3x2ά��Ԫ�������ÿ�б�ʾ����x��y����
	 * @param con
	 *            2x3ά���Գ���
	 * @param tp
	 * @param K
	 */
	public static void elementElasticityStiff(int[] ele, double[][] xy, double[][] con, BVPType tp, double[][] K) {
		double[] x = MVO.col(xy, 0);
		double[] N = new double[6];
		double[][] gradN = new double[6][2];
		double[] E = con[0];
		double[] Nu = con[1];
		double mE, mNu;
		double[] L = null;
		double[][] B = null;
		double area, weight;
		MVO.fill(K, 0.0);
		switch(tp) {
		case PLANESTRAIN:
		case PLANESTRESS:
			for(int i = 0; i < 4; i++) {
				L = Constants.TP3[i].clone();
				area = CG.T6ShapeFunction(xy, L, N, gradN);
				mE = MVO.dot_product(E, N);
				mNu = MVO.dot_product(Nu, N);
				double[][] D = Materials.ElasticConstituteMatrix(mE, mNu, tp);
				weight = area * Constants.TW3[i];
				B = new double[3][12];
				for(int ii = 0; ii < 6; ii++) {
					B[0][2 * ii] = B[2][2 * ii + 1] = gradN[ii][0];
					B[1][2 * ii + 1] = B[2][2 * ii] = gradN[ii][1];
				}
				double[][] KK = matmul(transpose(B), matmul(D, B));
				for(int si = 0; si < 12; si++) {
					for(int sj = 0; sj < 12; sj++) {
						K[si][sj] += weight * KK[si][sj];
					}
				}
			}
			break;
		case CONSTITUTE_MATRIX_COMMON:
			for(int i = 0; i < 3; i++) {
				L = Constants.TP2[i].clone();
				area = CG.T6ShapeFunction(xy, L, N, gradN);
				weight = area * Constants.TW2[i];
				B = new double[3][12];
				for(int ii = 0; ii < 6; ii++) {
					B[0][2 * ii] = B[2][2 * ii + 1] = gradN[ii][0];
					B[1][2 * ii + 1] = B[2][2 * ii] = gradN[ii][1];
				}
				double[][] KK = matmul(transpose(B), matmul(con, B));
				for(int si = 0; si < 12; si++) {
					for(int sj = 0; sj < 12; sj++) {
						K[si][sj] += weight * KK[si][sj];
					}
				}
			}
			break;
		case AXISYMMETRIC:
			for(int i = 0; i < 4; i++) {
				L = Constants.TP3[i].clone();
				area = CG.T6ShapeFunction(xy, L, N, gradN);
				mE = MVO.dot_product(E, N);
				mNu = MVO.dot_product(Nu, N);
				double[][] D = Materials.ElasticConstituteMatrix(mE, mNu, tp);
				weight = area * Constants.TW3[i] * MVO.dot_product(x, N);
				B = new double[4][12];
				for(int ii = 0; ii < 6; ii++) {
					B[0][2 * ii] = B[3][2 * ii + 1] = gradN[ii][0];
					B[1][2 * ii + 1] = B[3][2 * ii] = gradN[ii][1];
					B[2][2 * ii] = N[ii] / MVO.dot_product(x, N);
				}
				double[][] KK = matmul(transpose(B), matmul(D, B));
				for(int si = 0; si < 12; si++) {
					for(int sj = 0; sj < 12; sj++) {
						K[si][sj] += weight * KK[si][sj];
					}
				}
			}
			break;
		case CONSTITUTE_MATRIX_AXIS:
			for(int i = 0; i < 4; i++) {
				L = Constants.TP3[i].clone();
				area = CG.T6ShapeFunction(xy, L, N, gradN);
				mE = MVO.dot_product(E, N);
				mNu = MVO.dot_product(Nu, N);
				weight = area * Constants.TW3[i] * MVO.dot_product(x, N);
				B = new double[4][12];
				for(int ii = 0; ii < 6; ii++) {
					B[0][2 * ii] = B[3][2 * ii + 1] = gradN[ii][0];
					B[1][2 * ii + 1] = B[3][2 * ii] = gradN[ii][1];
					B[2][2 * ii] = N[ii] / MVO.dot_product(x, N);
				}
				double[][] KK = matmul(transpose(B), matmul(con, B));
				for(int si = 0; si < 12; si++) {
					for(int sj = 0; sj < 12; sj++) {
						K[si][sj] += weight * KK[si][sj];
					}
				}
			}
			break;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
	}

	public static void elementElasticityMass(int[] ele, double[][] xy, double[][] rho, BVPType tp, double[][] K) {
//		double[] p = rho[0];
//		double[][] grad = new double[3][2];
//		double area = CG.gradient2DT3(xy, grad);
//		double[] x = new double[] { xy[0][0], xy[1][0], xy[2][0] };
//		if (tp == BVPType.CONSTITUTE_MATRIX_AXIS || tp == BVPType.AXIS_LUMP_MATRIX) {
//			MVO.mul(p, x, p);
//		}
//		double N11 = (3 * p[0] + p[1] + p[2]) / 30 * area;
//		double N22 = (p[0] + 3 * p[1] + p[2]) / 30 * area;
//		double N33 = (p[0] + p[1] + 3 * p[2]) / 30 * area;
//		double N12 = (2 * p[0] + 2 * p[1] + p[2]) / 60 * area;
//		double N13 = (2 * p[0] + p[1] + 2 * p[2]) / 60 * area;
//		double N23 = (p[0] + 2 * p[1] + 2 * p[2]) / 60 * area;
//		MVO.fill(K, 0.0);
//		switch (tp) {
//		case CONSIST_MATRIX:
//		case AXIS_CONSIST_MATRIX:
//			K[0][0] = K[1][1] = N11;
//			K[0][2] = K[1][3] = N12;
//			K[0][4] = K[1][5] = N13;
//
//			K[2][0] = K[3][1] = N12;
//			K[2][2] = K[3][3] = N22;
//			K[2][4] = K[3][5] = N23;
//
//			K[4][0] = K[5][1] = N13;
//			K[4][2] = K[5][3] = N23;
//			K[4][4] = K[5][5] = N33;
//
//			break;
//		case LUMP_MATRIX:
//		case AXIS_LUMP_MATRIX:
//			K[0][0] = K[1][1] = (2 * p[0] + p[1] + p[2]) / 12 * area;
//			K[2][2] = K[3][3] = (p[0] + 2 * p[1] + p[2]) / 12 * area;
//			K[4][4] = K[5][5] = (p[0] + p[1] + 2 * p[2]) / 12 * area;
//			break;
//		default:
//			break;
//		}
		// System.out.println(MVO.toString(K));
		// System.exit(0);
	}

	public static void elementVolumeForce(int[] ele, double[][] xy, double[][] f, BVPType tp, double[] eleF) {
		double[] x = MVO.col(xy, 0);
		double[] N = new double[6];
		double[][] gradN = new double[6][2];
		double[] L = null;
		double mfu, mfv;
		double area, weight;
		Arrays.fill(eleF, 0.0);
		switch (tp) {
		case PLANESTRESS:
		case PLANESTRAIN:
		case CONSTITUTE_MATRIX_COMMON:
			for(int i = 0; i < 4; i++) {
				L = Constants.TP3[i].clone();
				area = CG.T6ShapeFunction(xy, L, N, gradN);
				mfu = MVO.dot_product(f[0], N);
				mfv = MVO.dot_product(f[1], N);
				weight = area  * Constants.TW3[i];
				for (int si = 0; si < 6; si++) {
					eleF[2 * si] += N[si] * weight * mfu;
					eleF[2 * si + 1] += N[si] * weight * mfv;
				}				
			}
			break;
		case AXISYMMETRIC:
		case CONSTITUTE_MATRIX_AXIS:
			for(int i = 0; i < 4; i++) {
				L = Constants.TP3[i].clone();
				area = CG.T6ShapeFunction(xy, L, N, gradN);
				mfu = MVO.dot_product(f[0], N);
				mfv = MVO.dot_product(f[1], N);
				weight = area  * Constants.TW3[i] * MVO.dot_product(x, N);
				for (int si = 0; si < 6; si++) {
					eleF[2 * si] += N[si] * weight * mfu;
					eleF[2 * si + 1] += N[si] * weight * mfv;
				}				
			}
			break;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
	}

	public static void elemenSurfaceForce(int[] ele, double[][] xy, double[][] f, BVPType tp, double[] eleF) {
		double length = CG.length2DL2(xy);
//		double[][] B = new double[6][6];
//		double[] fuv = {f[0][0], f[1][0], f[0][1], f[1][1], f[0][2], f[1][2]};
//		B[0][0] = B[1][1] = B[2][2] = B[3][3] = 2.0 * length / 15;
//		B[0][2] = B[2][0] = B[1][3] = B[3][1] = -length / 30.0;
//		B[0][4] = B[1][5] = B[2][4] = B[3][5] = length / 15.0;
//		B[4][0] = B[4][2] = B[5][1] = B[5][3] = length / 15.0;
//		B[4][4] = B[5][5] = 8.0 * length / 15.0;
//		MVO.matmul(B, fuv, ele);
		
		double[] gu = f[0];
		double[] gv = f[1];
		
		switch (tp) {
		case PLANESTRAIN:
		case PLANESTRESS:
		case CONSTITUTE_MATRIX_COMMON:
			eleF[0] = (4 * gu[0] - gu[1] + 2 * gu[2]) * length / 30;
			eleF[1] = (4 * gv[0] - gv[1] + 2 * gv[2]) * length / 30;
			eleF[2] = (-gu[0] + 4 * gu[1] + 2 * gu[2]) * length / 30;
			eleF[3] = (-gv[0] + 4 * gv[1] + 2 * gv[2]) * length / 30;
			eleF[4] = (gu[0] + gu[1] + 8 * gu[2]) * length / 15;
			eleF[5] = (gv[0] + gv[1] + 8 * gv[2]) * length / 15;
			break;
		case AXISYMMETRIC:
		case CONSTITUTE_MATRIX_AXIS:
			double[] x = MVO.col(xy, 0);
			eleF[0] = (4 * gu[0] * x[0] - gu[1] * x[1] + 2 * gu[2] * x[2]) * length / 30;
			eleF[1] = (4 * gv[0] * x[0] - gv[1] * x[1] + 2 * gv[2] * x[2]) * length / 30;
			eleF[2] = (-gu[0] * x[0] + 4 * gu[1] * x[1] + 2 * gu[2] * x[2]) * length / 30;
			eleF[3] = (-gv[0] * x[0] + 4 * gv[1] * x[1] + 2 * gv[2] * x[2]) * length / 30;
			eleF[4] = (gu[0] * x[0] + gu[1] * x[1] + 8 * gu[2] * x[2]) * length / 15;
			eleF[5] = (gv[0] * x[0] + gv[1] * x[1] + 8 * gv[2] * x[2]) * length / 15;
			break;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
//		double[] gu = f[0];
//		double[] gv = f[1];
//		double sumGU = MVO.sum(gu);
//		double sumGV = MVO.sum(gv);
//		switch (tp) {
//		case PLANESTRAIN:
//		case PLANESTRESS:
//		case CONSTITUTE_MATRIX_COMMON:
//			length /= 6.0;
//			for (int i = 0; i < 2; i++) {
//				ele[2 * i] = (sumGU + gu[i]) * length;
//				ele[2 * i + 1] = (sumGV + gv[i]) * length;
//			}
//			break;
//		case AXISYMMETRIC:
//		case CONSTITUTE_MATRIX_AXIS:
//			length /= 12.0;
//			double[] x = MVO.col(xy, 0);
//			double sumX = MVO.sum(x);
//			for (int i = 0; i < 2; i++) {
//				ele[2 * i] = (sumGU * sumX + 2 * gu[i] * x[i]) * length;
//				ele[2 * i + 1] = (sumGV * sumX + 2 * gv[i] * x[i]) * length;
//			}
//			break;
//		default:
//			throw new IllegalArgumentException("Wrong BVPType!");
//		}
	}

	public void assemblePressure(ScalarFunc[] coefFunc, double[] param, BVPType tp, double[] RHS, int[] bdLabel) {
		assembleBoundaryVector(coefFunc, param, FES2T62::elementPressure, tp, RHS, bdLabel);
	}

	public void assemblePressure(double[] constCoef, BVPType tp, double[] RHS, int[] bdLabel) {
		assembleBoundaryVector(constCoef, FES2T62::elementPressure, tp, RHS, bdLabel);
	}
	
	public void assemblePressure(double[][] vecCoef, BVPType tp, double[] RHS, int[] bdLabel) {
		assembleBoundaryVector(vecCoef, FES2T62::elementPressure, tp, RHS, bdLabel);
	}
	
	public static void elementPressure(int[] ele, double[][] xy, double[][] f, BVPType tp, double[] eleF) {
		double[] pressure = f[0];
		double[] x = MVO.col(xy, 0);
		double[] y = MVO.col(xy, 1);
		double len1 = Math.sqrt((x[0] - x[2])*(x[0] - x[2]) + (y[0] - y[2]) * (y[0] - y[2]));
		double nx1 = (y[2] - y[0])/len1;
		double ny1 = (x[0] - x[2])/len1;
		double len2 = Math.sqrt((x[1] - x[2])*(x[1] - x[2]) + (y[1] - y[2]) * (y[1] - y[2]));
		double nx2 = (y[1] - y[2])/len2;
		double ny2 = (x[2] - x[1])/len2;
		f[0][0] = -pressure[0] * nx1;
		f[0][1] = -pressure[1] * nx1;
		f[0][2] = -pressure[2] * 0.5 * (nx1 + nx2);
		f[1][0] = -pressure[0] * ny2;
		f[1][1] = -pressure[1] * ny2;
		f[1][2] = -pressure[2] * 0.5 * (ny1 + ny2);
		elemenSurfaceForce(ele, xy, f, tp, eleF);
	}	

	public void assembleHeatLoad(ScalarFunc[] coef, double[] param, double[][] vecCoef, BVPType tp, double[] RHS) {
		assembleGlobalVector(coef, param, vecCoef, FES2T62::elementHeatLoad, tp, RHS);
	}

	public void assembleHeatLoad(double[] coef, double[][] vecCoef, BVPType tp, double[] RHS) {
		assembleGlobalVector(coef, vecCoef, FES2T62::elementHeatLoad, tp, RHS);
	}

	public static void elementHeatLoad(int[] ele, double[][] coord, double[][] coef, BVPType tp, double[] eleLoad) {
//		double[] beta = coef[0];
//		double[] temperature = coef[1];
//		double[] x = MVO.col(coord, 0);
//		double[][] grad = new double[3][2];
//		double area = CG.gradient2DT3(coord, grad);
//		Arrays.fill(eleLoad, 0.0);
//		switch (tp) {
//		case PLANESTRAIN:
//		case PLANESTRESS:
//		case CONSTITUTE_MATRIX_COMMON:
//			double c = (MVO.sum(beta) * MVO.sum(temperature) + MVO.dot_product(beta, temperature)) * area / 12.0;
//			for (int i = 0; i < 3; i++) {
//				eleLoad[2 * i] = c * grad[i][0];
//				eleLoad[2 * i + 1] = c * grad[i][1];
//			}
//			break;
//		case AXISYMMETRIC:
//		case CONSTITUTE_MATRIX_AXIS:
//			for (int i = 0; i < 4; i++) {
//				double[] L = Constants.TP3[i];
//				double wxarea = MVO.dot_product(beta, L) * MVO.dot_product(temperature, L) * TW3[i];
//				for (int j = 0; j < 3; j++) {
//					eleLoad[2 * j] += (grad[j][0] * MVO.dot_product(x, L) + L[j]) * wxarea;
//					eleLoad[2 * j + 1] += grad[j][1] * MVO.dot_product(x, L) * area;
//				}
//			}
//			break;
//		default:
//			throw new IllegalArgumentException("Wrong BVPType!");
//		}
		// System.out.println(Arrays.toString(eleLoad));

	}

	public void computeNodalStrain(double[] sol, double[][] strain) {
//		int len = sol.length;
//		ParamCheck.checkEqual(len, this.ndof);
//		int n = len / 2;
//		ParamCheck.checkEqual(n, strain[0].length);
//		ParamCheck.checkEqual(n, strain[1].length);
//		ParamCheck.checkEqual(n, strain[2].length);
//		double[] angle = new double[3];
//		double[][] coord = new double[3][2];
//		double[][] grad = new double[3][2];
//		double[] angleCount = new double[this.GetNdof()];
//		double dxU, dyU, dxV, dyV;
//		Element ele = null;
//		double[] eleU = new double[3];
//		double[] eleV = new double[3];
//		for (int i = 0; i < mesh.getNt(); i++) {
//			ele = mesh.getElement(i);
//			mesh.getCoordPerElement(ele, coord);
//			CG.gradient2DT3(coord, grad);
//			CG.angle2DT3(coord, angle);
//			for (int j = 0; j < 3; j++) {
//				eleU[j] = sol[2 * ele.idx[j]];
//				eleV[j] = sol[2 * ele.idx[j] + 1];
//			}
//			double[] dx = MVO.col(grad, 0);
//			double[] dy = MVO.col(grad, 1);
//			dxU = MVO.dot_product(eleU, dx);
//			dyU = MVO.dot_product(eleU, dy);
//			dxV = MVO.dot_product(eleV, dx);
//			dyV = MVO.dot_product(eleV, dy);
//			for (int j = 0; j < 3; j++) {
//				strain[0][ele.idx[j]] += dxU * angle[j];
//				strain[1][ele.idx[j]] += dyV * angle[j];
//				strain[2][ele.idx[j]] += (dxV + dyU) * angle[j];
//				angleCount[ele.idx[j]] += angle[j];
//			}
//		}
//		for (int i = 0; i < n; i++) {
//			strain[0][i] /= angleCount[i];
//			strain[1][i] /= angleCount[i];
//			strain[2][i] /= angleCount[i];
//		}
		
//		int len = u.length;
//		ParamCheck.checkEqual(len, this.ndof);
//
//		double[] angle = new double[6];
//		double[][] coord = new double[6][2];
//		double[][] grad = new double[3][2];
//		double[][] eleGradU = new double[6][2];
//		double[] eleSDU = new double[4];
//		double[] angleCount = new double[this.GetNdof()];
//		Element ele = null;
//		double[] eleU = new double[6];
//		for(int i = 0, n = mesh.getNt(); i < n; i++) {
//			ele = mesh.getElement(i);
//			mesh.getCoordPerElement(ele, coord);
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = u[ele.idx[j]];
//			}
//			CG.gradient2DT3(coord, grad);
//			CG.gradientValue2DT6(grad, eleU, eleGradU);
//			CG.secondDerivative2DT6(grad, eleU, eleSDU);
//			CG.angle2DT6(coord, angle);
//			for(int j = 0; j < 6; j++) {
//				fdU[0][ele.idx[j]] += eleGradU[j][0] * angle[j];
//				fdU[1][ele.idx[j]] += eleGradU[j][1] * angle[j];
//				sdU[0][ele.idx[j]] += eleSDU[0] * angle[j];
//				sdU[1][ele.idx[j]] += eleSDU[1] * angle[j];
//				sdU[2][ele.idx[j]] += eleSDU[2] * angle[j];
//				angleCount[ele.idx[j]] += angle[j];
//			}
//		}
//		for(int i = 0, n = mesh.getNv(); i < n; i++) {
//			fdU[0][i] /= angleCount[i];
//			fdU[1][i] /= angleCount[i];
//			sdU[0][i] /= angleCount[i];
//			sdU[1][i] /= angleCount[i];
//			sdU[2][i] /= angleCount[i];
//		}
//		
//		for(int i = 0; i < mesh.getNt(); i++) {
//			ele = mesh.getElement(i);
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = sdU[0][ele.idx[j]];
//			}
//			sdU[0][ele.idx[3]] = 0.5 * (eleU[0] + eleU[1]);
//			sdU[0][ele.idx[4]] = 0.5 * (eleU[1] + eleU[2]);
//			sdU[0][ele.idx[5]] = 0.5 * (eleU[0] + eleU[2]);
//			
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = sdU[1][ele.idx[j]];
//			}
//			sdU[1][ele.idx[3]] = 0.5 * (eleU[0] + eleU[1]);
//			sdU[1][ele.idx[4]] = 0.5 * (eleU[1] + eleU[2]);
//			sdU[1][ele.idx[5]] = 0.5 * (eleU[0] + eleU[2]);
//			
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = sdU[2][ele.idx[j]];
//			}
//			sdU[2][ele.idx[3]] = 0.5 * (eleU[0] + eleU[1]);
//			sdU[2][ele.idx[4]] = 0.5 * (eleU[1] + eleU[2]);
//			sdU[2][ele.idx[5]] = 0.5 * (eleU[0] + eleU[2]);
//		}
	}


	/**
	 * 
	 * @param uv
	 * @param fdUV
	 * @param sdUV
	 */
	public void computeFirstAndSecondDerivative(double[][] uv, double[][] fdUV, double[][] sdUV) {
		int len = uv[0].length;
		ParamCheck.checkEqual(2 * len, this.ndof);

		double[] angle = new double[6];
		double[][] coord = new double[6][2];
		double[][] grad = new double[3][2];
		double[][] eleGradU = new double[6][2];
		double[][] eleGradV = new double[6][2];
		double[] eleSDU = new double[3];
		double[] eleSDV = new double[3];
		double[] angleCount = new double[this.GetNdof()];
		double[] eleU = new double[6];
		double[] eleV = new double[6];
		for(int i = 0, n = mesh.getNt(); i < n; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coord);
			for(int j = 0; j < 6; j++) {
				eleU[j] = uv[0][ele[j]];
				eleV[j] = uv[1][ele[j]];
			}
			CG.gradient2DT3(coord, grad);
			CG.gradientValue2DT6(grad, eleU, eleGradU);
			CG.gradientValue2DT6(grad, eleV, eleGradV);
			CG.secondDerivative2DT6(grad, eleU, eleSDU);
			CG.secondDerivative2DT6(grad, eleV, eleSDV);
			CG.angle2DT6(coord, angle);
			for(int j = 0; j < 6; j++) {
				fdUV[0][ele[j]] += eleGradU[j][0] * angle[j];
				fdUV[1][ele[j]] += eleGradU[j][1] * angle[j];
				fdUV[2][ele[j]] += eleGradV[j][0] * angle[j];
				fdUV[3][ele[j]] += eleGradV[j][1] * angle[j];
				sdUV[0][ele[j]] += eleSDU[0] * angle[j];
				sdUV[1][ele[j]] += eleSDU[1] * angle[j];
				sdUV[2][ele[j]] += eleSDU[2] * angle[j];
				sdUV[3][ele[j]] += eleSDV[0] * angle[j];
				sdUV[4][ele[j]] += eleSDV[1] * angle[j];
				sdUV[5][ele[j]] += eleSDV[2] * angle[j];
				angleCount[ele[j]] += angle[j];
			}
		}
		for(int i = 0, n = mesh.getNv(); i < n; i++) {
			fdUV[0][i] /= angleCount[i];
			fdUV[1][i] /= angleCount[i];
			fdUV[2][i] /= angleCount[i];
			fdUV[3][i] /= angleCount[i];
			sdUV[0][i] /= angleCount[i];
			sdUV[1][i] /= angleCount[i];
			sdUV[2][i] /= angleCount[i];
			sdUV[3][i] /= angleCount[i];
			sdUV[4][i] /= angleCount[i];
			sdUV[5][i] /= angleCount[i];
		}
		
		for(int i = 0; i < mesh.getNt(); i++) {
			int[] ele = mesh.getElement(i);
			for(int k = 0; k < 6; k++) {
				for(int j = 0; j < 6; j++) {
					eleU[j] = sdUV[k][ele[j]];
				}
				sdUV[k][ele[3]] = 0.5 * (eleU[0] + eleU[1]);
				sdUV[k][ele[4]] = 0.5 * (eleU[1] + eleU[2]);
				sdUV[k][ele[5]] = 0.5 * (eleU[0] + eleU[2]);
			}
		}
			
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = sdUV[1][ele.idx[j]];
//			}
//			sdUV[1][ele.idx[3]] = 0.5 * (eleU[0] + eleU[1]);
//			sdUV[1][ele.idx[4]] = 0.5 * (eleU[1] + eleU[2]);
//			sdUV[1][ele.idx[5]] = 0.5 * (eleU[0] + eleU[2]);
//			
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = sdUV[2][ele.idx[j]];
//			}
//			sdUV[2][ele.idx[3]] = 0.5 * (eleU[0] + eleU[1]);
//			sdUV[2][ele.idx[4]] = 0.5 * (eleU[1] + eleU[2]);
//			sdUV[2][ele.idx[5]] = 0.5 * (eleU[0] + eleU[2]);
//
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = sdUV[0][ele.idx[j]];
//			}
//			sdUV[0][ele.idx[3]] = 0.5 * (eleU[0] + eleU[1]);
//			sdUV[0][ele.idx[4]] = 0.5 * (eleU[1] + eleU[2]);
//			sdUV[0][ele.idx[5]] = 0.5 * (eleU[0] + eleU[2]);
//			
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = sdUV[1][ele.idx[j]];
//			}
//			sdUV[1][ele.idx[3]] = 0.5 * (eleU[0] + eleU[1]);
//			sdUV[1][ele.idx[4]] = 0.5 * (eleU[1] + eleU[2]);
//			sdUV[1][ele.idx[5]] = 0.5 * (eleU[0] + eleU[2]);
//			
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = sdUV[2][ele.idx[j]];
//			}
//			sdUV[2][ele.idx[3]] = 0.5 * (eleU[0] + eleU[1]);
//			sdUV[2][ele.idx[4]] = 0.5 * (eleU[1] + eleU[2]);
//			sdUV[2][ele.idx[5]] = 0.5 * (eleU[0] + eleU[2]);
		
	}
	
}
